#include <Element.d/Spring.d/LinSpring.h>
#include <Utils.d/dofset.h>
#include <Utils.d/linkfc.h>
#include <Corotational.d/SpringCorotator.h>
#include <Math.d/FullSquareMatrix.h>

extern "C"      {
void    _FORTRAN(mstf11)(double*,double&,double&,double&,double&,double&,
                         double&,double&,double&,double&,double&,double&,
                         double&);
}

LinSpring::LinSpring(int* nodenums)
{
        nn[0] = nodenums[0];
}

Element *
LinSpring::clone()
{
 return new LinSpring(*this);
}

void
LinSpring::renum(int *table)
{
  nn[0] = table[nn[0]];
}

void
LinSpring::renum(EleRenumMap& table)
{
  nn[0] = table[nn[0]];
}

FullSquareMatrix
LinSpring::massMatrix(CoordSet &,double *mel,int cmflg)
{
        FullSquareMatrix elementMassMatrix(3,mel);

	// zero the element mass matrix
	elementMassMatrix.zero();

        return elementMassMatrix;
}

FullSquareMatrix
LinSpring::stiffness(CoordSet &, double *d, int flg)
{
        double sc1 = prop->A;
        double sc2 = prop->E;
        double sc3 = prop->nu;

        double u1  = prop->rho;
        double u2  = prop->c;
        double u3  = prop->k;

        double v1  = prop->eh;
        double v2  = prop->P;
        double v3  = prop->Ta;

        double w1  = prop->Q;
        double w2  = prop->W;
        double w3  = prop->Ixx;

	// KHP Sometime modify this so that there is no call to FORTRAN
	// the stiffness matrix is easy to do right here in C++

	_FORTRAN(mstf11)((double*)d,sc1,sc2,sc3,u1,u2,u3,v1,v2,v3,w1,w2,w3);

        FullSquareMatrix ret(3,d);

        return ret;
}

int
LinSpring::numNodes()
{
        return 1;
}

int *
LinSpring::nodes(int *p)
{
        if(p == 0) p = new int[1];
        p[0] = nn[0];
        return p;
}

int
LinSpring::numDofs()
{
        return 3;
}

int *
LinSpring::dofs(DofSetArray &dsa, int *p)
{
        if(p == 0) p = new int[3];
      
	dsa.number(nn[0],DofSet::XYZdisp,p);

        return p;
}

void
LinSpring::markDofs( DofSetArray &dsa )
{
        dsa.mark( nn[0], DofSet::XYZdisp );
}

Corotator *
LinSpring::getCorotator(CoordSet &cs, double *kel, int , int )
{
 int flag = 0;
 int numnodes = 1;
 FullSquareMatrix myStiff = stiffness(cs, kel, flag);

 return new SpringCorotator( nn, cs, numnodes, prop->A, prop->E, prop->nu,
                                myStiff);

}

